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Abstract — Qiiasi-onedimensional stereoregular polymers as for example polyacetylene 
are currently of considerable interest, both experimentally as well as theoretically. There are 
basically two different approaches for doing electronic structure calculations: One method, 
the so-called crystal orbital method, uses periodic boundary conditions and is essentially 
based on concepts of solid state theory. The other method is essentially a quantum chemical 
method since it approximates the polymer by oligomers consisting of a finite number of 
monomer units, i.e., by molecules of finite size. In this way, the highly developed technology 
of quantum chemical molecular programs can be used. Unfortunately, oligomers of finite 
size are not necessarily able to model those features of a polymer which crucially depend 
of its in principle infinite extension. However, in such a case extrapolation techniques can 
be extremely helpful. For example, one can perform electronic structure calculations for a 
sequence of oligomers with an increasing number of monomer units. In the next step, one 
then can try to determine the limit of this sequence for an oligomer of infinite length with 
the help of suitable extrapolation methods. Several different extrapolation methods are dis- 
cussed which are able to accomplish an extrapolation of energies and properties of oligomers 
to the infinite chain limit. Calculations for the ground state energy of polyacetylene are 
presented which demonstrate the practical usefulness of extrapolation methods. 

Keywords — Quantum chemical oligomer calculations, quasi-onedimensional stereoregu- 
lar polymers, extrapolation methods, sequence transformations. 



1 Introduction 

The classic example of a gwsi-onedimensional stereoregular polymer is polyacetylene, which 
can be generated by translating the repeat unit — (CH=CH) — . Conjugated polymers such as 
polyacetylene are currently of great technological interest because of potential applications in 
a wide variety of optical and optoelectronic devices based on their conducting, photophysical, 
and nonlinear optical properties Moreover, the keen interest in polyacetylene especially is 

documented by the fact that there even exists a monograph with the title Polyacetylene [0] . 

Numerous properties involving the electronic structure of these polymers are not yet under- 
stood in a completely satisfactory way. Consequently, these polymers are not only very interesting 
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from a technological point of view but are also quite challenging from a basic science point of 
view. 

There are two approaches to ab initio molecular electronic structure calculations. One may 
be described as the wavefunction approach and the other as the density approach (also known 
as density functional theory). The starting point for the wavefunction approach is an effective 
independent particle treatment of the Schrodinger equation named after Hartree Q and Fock H . 
In the density approach one uses an equation of similar general form due to Kohn and Sham [|6| 
but, in principle, the electron interactions are fully taken into account. The solution of cither 
equation is a set of one-electron functions, called orbitals, that can be used to construct the multi- 
electron wavefunction (Hartree- Fock) or the electron density (Kohn- Sham). In cither case the 
orbitals are expanded in terms of finite linear combination of analytic basis functions (nowadays 
usually Gaussians). 

The Hartree-Fock approach combined with analytic basis functions leads to the so-called 
Roothaan equations which can be written in the form of a generalized matrix eigenvalue 
problem. For meaningful results it is necessary for the number of basis functions and, hence, the 
dimension of the Roothaan equations to exceed the number of electrons. This imposes obvious 
restrictions on the size of the size of oligomers used to model a polymer, which can be handled 
by ab initio programs based on the highly developed computational techniques of molecular 
quantum chemistry (see for instance ||,||]). 

If, however, the polymer under consideration is a stereoregular qwsi-onedimensional polymer 
like polyacetylene, it is nevertheless possible to compute energies or properties on the basis of 
quantum mechanics. There are two fundamentally different computational approaches - the 
crystal orbital approach and the finite oligomer method - that have complementary advantages 
and disadvantages. 

The crystal orbital approach, which is for example discussed in monographs by Pisani, Dovesi, 
and Roetti p0| , Ladik (TTJ , and Andre, Delhalle, and Bredas or in review articles by Andre 
p3[ , Kertesz Q , and Ladik [|l5| , treats polymers as solids and uses the concepts and techniques 
of solid state physics such as translational symmetry and cyclic boundary conditions (see for 
example [^6[0). This leads to band structure, or crystal orbital, calculations whose complexity is 
limited by the number of interacting neighboring unit cells. Such calculations have the undeniable 
advantage that one only has to consider these interactions for a representative unit cell since all 
unit cells are identical. This is a major simplification over a molecular-type calculation where 
the unit cells are translationally inequivalent. Thus, the crystal orbital approach is in principle 
well suited for describing those features of a polymer which crucially depend on its in principle 
infinite extension. 

Nevertheless, there are still some open problems with crystal orbital calculations. For ex- 
ample, convergence problems in the Hartree-Fock iterations and the proper handling of long 
range Coulomb interactions can be problematic jL8|. Moreover, it is very difficult to take into 
account the effect of nonperiodic perturbations such as the interaction with an electric field. The 
term in the polarization response, which is linear in the interacting electric field, determines the 
polarizability of the polymer, whereas the nonlinear terms determine the hyperpolarizabilities. 
The nonlinear response is connected with properties such as second harmonic generation, optical 
rectification, and the intensity dependent refractive index (see for example [^9[) which are of 
considerable interest in materials science these days. Despite many years of effort going back 
to the 1940's and substantial recent progress p0|-p5[, it is fair to say that an incontrovertible 
crystal orbital method |2(| has yet to be presented, even at a simplified Hartree-Fock level 
where just the average - rather than the instantaneous - interelectronic repulsion is taken into 
account. 
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Another problem of the crystal orbital approach is that the explicit incorporation of electronic 
correlation is difficult. In the Hartree-Fock approximation the instantaneous interaction between 
electrons is replaced by an orbital-averaged interaction. The error due to making this approxi- 
mation is known as the electron correlation error (sec for example [ 28 - 30 1 ) . For some properties, 
such as the energy per unit cell, this error is quite small on a percentage basis. Nonetheless, 
correlation is very important because it can easily be the deciding factor between alternative 
structural conformations. For other properties, such as hypcrpolarizabilities, the correlation 
error may be as large as 800 % 

The density method includes electron correlation but the one-electron potential, which simul- 
taneously takes into account both electron exchange (Pauli Principle) and correlation, must be 
approximated. In practice, the approximations used are local in nature - they involve the den- 
sity, its gradient and, possibly, the Laplacian at each point. For that reason the method fails 
to predict electrical properties (dipole moment, linear polarizability, hyperpolarizabilites) ]3^,^| 
of spatially extended systems. One can expect that the same will be true of other properties 
of spatially extended systems, such as the mechanical force constants describing the interaction 
between different unit cells, but that has not been checked as yet. 

Because of the problems mentioned above, a different approach - the finite oligomer method 
- has gained popularity in the middle to the late 1980's: Stereoregular guasi-onedimensional 
polymers are approximated by molecules of finite length []35|-|42f . In the case of polyacetylene, this 
means that the polymer (013=011)00 is approximated by an oligomer H— (CH=CH)n— H, where 
TV is a finite integer. The oligomer approach offers many advantages. The ab initio treatment of 
electron correlation using any of the techniques available in standard molecular ab initio programs 
is straightforward. Moreover, the same is true for the calculation of hyperpolarizabilities, which 
has yet to be successfully implemented within the framework of crystal orbital theory even at 
the Hartree-Fock level. Molecular-type calculations using finite oligomers have been done not 
only at the Hartree-Fock level, but have also incorporated electron correlation by many-body 
perturbation theory S-B. As a result, it has become known that the role of electron 
correlation can be exceptionally significant for this problem. The same may be said regarding 
the importance of the vibrational contribution to the hyperpolarizability [^6|, ^| that arises 
because the electronic response depends upon the nuclear positions. In addition, studies have 
been undertaken to determine the effect of interactions in the solid state with other polymer 
chains |^8|, fl9| and with dopants [j50|~j52| that either accept or donate electrons. All of these 
investigations are readily feasible because of the availability of quantum chemistry computer 
codes developed for ordinary finite molecules. 

However, the finite oligomer approach of course also has some shortcomings which may lead 
to numerical problems. The basic problem of the approximation of a quasi-lineax stereoregular 
polymer Aoo by an oligomer X — An— Y, whose terminal free valences are saturated by suitably 
terminal groups X and Y (usually hydrogen atoms), is the choice of the number N of repeat units 
A. Of course, N should be small enough to permit an accurate quantum chemical calculation of 
the oligomer. However, N should also be large enough in order to provide a sufficiently realistic 
representation of those features of the polymer that result from its in principle infinite extension. 

Obviously, these two requirements tend to be contradictory which implies that finite oligomer 
calculations can suffer badly from slow convergence to the infinite chain limit. This is where 
extrapolation methods come into play: Instead of calculating the intensive energy E(N) or 
more generally for an intensive property P{N) of a single oligomer X— An— Y, one can perform 
calculations for a sequence of oligomers X— An— Y with N = 1, 2, . . . , N max . The sequences E(l), 
E(2), . . . , E(N mayi ) and P(l), -P(2), . . . , P(N max ) can then be extrapolated to the infinite chain 
limits E(oo) and P(oo), respectively. 
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Many different approaches for doing the extrapolations are possible. In this article, we want 
to use sequence transformations for that purpose. Numerical techniques for the acceleration of 
convergence or the summation of divergent series are as old as calculus. According to Knopp fl53| , 
p. 249], the first series transformation was published by Stirling (54j already in 1730, and in 1755 
Euler published the series transformation which now bears his name. In rudimentary form, 
convergence acceleration methods are even older. Brezinski |p6| pp. 90 - 91] mentioned that in 
1654 Huygens used a linear extrapolation scheme which is a special case of what we now call 
Richardson extrapolation ]57j for obtaining better approximations to 7r, and that in 1674 Seki 
Kowa, the probably most famous Japanese mathematician of that period, tried to accomplish 
this with the help of the so-called A 2 process, which is usually attributed to Aitken J58|. Further 
details on the historical development of extrapolation methods can be found in an article by 
Brezinski @. 

The modern theory of sequence transformations begins with two articles by Shanks |60) and 
Wynn jtnjl . Shanks introduced a sequence transformation which in the case of a power series 
produces Pade approximants |6^ | , and Wynn showed that this transformation can be computed 
conveniently by a nonlinear recursive scheme called the epsilon algorithm. These two articles 
appeared at a time when the first computers became publicly available. The subsequent revo- 
lution created a great demand for powerful numerical tools, and the two articles by Shanks and 
Wynn, which offered solutions to practical problems, had an enormous impact: They stimulated 
research not only on Pade approximants ]62|] but also on other sequence transformations, as doc- 
umented by the large number of monographs p3|-[72|] and review articles [|73|,[74| which appeared 
in the last years twenty years or so. Moreover, there is not only a lot of research on sequence 
transformations and their properties, but the number of articles which describe applications also 
increases steadily (see for example the list of references in j75| ) . 

Sequence transformations have already been used for the extrapolation of oligomer calcula- 
tions. Initially, they were applied to ground state energies but more recently, they have 
also been employed for the the computation of static and dynamic linear polarizabilities ]76[ ] . It 
is the intention of this article to provide a sufficiently detailed discussion of those sequence trans- 
formations that are suited for the extrapolation of oligomer calculations to the infinite limits. 
As discussed later in more detail, sequence transformations try to accomplish an improvement 
of convergence by detecting and utilizing regularities of the input data. For that purpose, they 
have to access the information stored in the later digits of the input data. However, in the case 
of oligomer calculations, the input data are produced by molecular ab initio programs which are 
huge packages of FORTRAN code (more than 100 000 lines of code). Normally, these programs 
operate in FORTRAN DOUBLE PRECISION which corresponds to an accuracy of 14 - 16 dec- 
imal digits (depending on the compiler). Unfortunately, this does not mean that their results 
have this accuracy. Numerous internal approximations plus the inevitable rounding errors reduce 
the accuracy of the results. This alone is not necessarily a problem. However, the complexity 
of the calculations, which are performed in such a molecular program, makes it impossible to 
obtain an estimate of the errors by standard mathematical approaches, and the only realistic 
alternative are mathematical experiments. Consequently, the emphasis of this manuscript will 
not be so much on powerful sequence transformations but rather on sequence transformations 
which produce reliable results under the restrictions that are typical of oligomer calculations. 

In Section (|^) , we give a short discussion of use the Euler Maclaurin formula for the evaluation 
of infinite series and show that sequence transformations also try to eliminate the truncations 
errors, albeit by purely numerical means. In Sections (|) and (|), we discuss sequence transfor- 
mations that are able to acelerate either linear or logarithmic convergence. In Section (^|), we 
show that oligomer calculations for the grond state energy of polyacetylene can be extrapolated 
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effectively to the infinite chain limit by sequence transformations, and finally in Section (|6|) there 
is a short summary. 



2 The Euler-Maclaurin Formula and Sequence Transfor- 
mations 

A model example of a slowly convergent series is the Dirichlet series for the Riemann zeta 
function: 

oo 

C(z) = 5> + l)-*. (2.1) 

This series converges provided that Re(z) > 1. However, it is notorious for extremely slow 
convergence if Re(z) is only slightly larger than 1. 

The Euler-Maclaurin formula, which is treated in most books on asymptotics (see for instance 
\ \n[ pp. 279 - 295]), makes it possible to evaluate C(z) even if the convergence of the Dirichlet 
series is extremely slow. Let us consider an infinite series X^l^Lo ana - ^ us assume that its 
terms f{v) are smooth and slowly varying functions of the index v. Then, the integral 

N 

f(x)dx (2.2) 

M 

with M, N G Z can provide a good approximation to the sum 

i/(Af) + f(M + 1) + . . . + f(N - 1) + i/(A0 (2.3) 

and vice versa. To see this, one only has to approximate the integrand by step functions. 

In the years between 1730 and 1740, Euler and Maclaurin derived independently correction 
terms, which ultimately yielded what we now call the Euler-Maclaurin formula: 



£ f[y) = / f(x)dx + [f(M) + f(N)} 

,=M Jm 1 



1 Bo 



E (gy? [f^HN) fto-VMl + R k (f) , (2.4a) 



Rk{}) = ~j2ky. j M B2k{ - X ~ W)/ (2fe) Wdx. (2.4b) 



Here, \x\ is the integral part of x, Bk(x) is a Bernoulli polynomial , and B m = B rn (0) is a 
Bernoulli number. 

If we set M = n + 1 and N = oo, the leading terms of the Euler-Maclaurin formula provides 
rapidly convergent approximations to the truncation error X)^Ln+i/(^)i provided that / is 
smooth and a slowly varying function of the index v. Of course, it is also necessary that the 
integral as well as derivatives of / with respect to v can be computed. 

The terms [y + 1)~ 2 of the Dirichlet series ( pTl| ) are obviously smooth functions of the index 
v and they can be differentiated and integrated easily. Thus, we can apply the Euler-Maclaurin 
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formula (2.4) with M = n+ 1 and N = oo to the truncation error of the Dirichlct scries: 



, s , (n + 2) 1 - z 1 , 



z - 1 

+ E (» + 2 )"* _SM+1 + ^ 2 ) ' ( 2 - 5a ) 



(£)» r B 2k ( x - ix}) 

(2k)\ J n+1 (1 + xY 



R k(n,z) = / n , ^ z+ 2k dx - ( 2 - 5b ) 



Here, (z) m — z(z + 1) • • • (z + m — 1) = T(z + m)/r(z) is a Pochhammer symbol. 



The expansion (2.5) for the truncation error of the Dirichlct scries (2.1) for £(z) is only 
asymptotic as n — > oo. Nevertheless, it is numerically extremely useful for sufficiently large 
values of n, as shown by the following simple example: 

C(«) - E(^ + ir + + \ (n + 2)- + E ^^^^ (» + 2)- 2 " 2j - +1 ■ (2-6) 

If we now choose n — 20 and z = 1.01, we obtain 

n 

^T(v + 1)- Z = 3.599 497 439 829 47, (2.7a) 

96.956 241 819 2202, (2.7b) 
0.220 355 095 043 682 x 10" 1 , (2.7c) 
0.168 605 034 844 029 x 10~ 3 , (2.7d) 



u=0 

(n 



- 1 
1 



2(n + 2) 2 
(z)iB 2 



2\(n + 2) z+1 

' ' ' :: L 1 - -0.351 266 295 216 895 x 10~ 7 , (2.7e) 



4!(n + 2) 2 + 3 



!: ' L "' - 0.347 155 401 295 600 x 10~ 10 . (2.7f) 



10!(n + 2) 2 + 5 



The sum of these 6 terms, which is 100.577 943 338 497, agrees completely with the "exact" 
result C(l-01) = 100.577 943 338 497 obtained by the computer algebra system Maple. 



This numerical example shows that a suitable truncation of the Euler-Maclaurin formula (2.4) 
provides an excellent approximation to the truncation error X)^n+i( t/ + -0 ~ z ■ Here, it should be 
taken into account that it is practically impossible to compute £(1.01) with sufficient accuracy 
by adding up the terms of the Dirichlet series (2.1) successively. We would need n « 10 600 to 



obtain (n + 2) 1 ^ z j[z — 1) = 10~ 4 for the leading term on the right-hand side of Eq. ( |2.4|) , which 
corresponds to the integral ([2.2]). Thus, 10 600 terms of the Dirichlet series ( |2.l| ) with z = 1.01 
would lead to an accuracy of only 6 decimal places in the final result. 

In view of these excellent results, it looks like an obvious idea to try to use the Euler-Maclaurin 
formula also in the case of other slowly convergent series. Unfortunately, this is not always 
possible. The Euler-Maclaurin formula requires that the terms of the series can be differentiated 
and integrated with respect to the index. This is only possible if the terms have a sufficiently 
simple analytic structure, which excludes most cases of interest. Moreover, the Euler-Maclaurin 
formula cannot be applied in the case of alternating or divergent series since their terms are 
neither smooth nor slowly varying functions of the index. 

However, the probably worst drawback of the Euler-Maclaurin formula is that it is an analytic 
convergence acceleration method. This means that it cannot be applied if only the numerical 



Extrapolation Methods for Oligomer Calculations 



7 



values of the terms of a series are known. This is quite unfortunate since convergence acceleration 
techniques are often badly needed in situations in which apart from a few numerical data very 
little is known. 

These facts give rise to the obvious question whether the Euler-Maclaurin formula can be 
generalized in such a way that approximations to the truncations errors of sequences and series 
can be constructed under less restrictive conditions. 

Principal tools to overcome convergence problems are sequence transformations which can 
accomplish something even if no explicit analytical expression for the truncation error is available. 
Here, it should be noted that Pade approximants, which have become the standard tool to 
overcome convergence problems with power series, are special sequence transformations since 
the partial sums of a power series are transformed into a doubly indexed sequence of rational 
functions J62[ . 

Let us assume that we have a sequence {s„}^ =0 which converges to some limit s = Soc. Then, 
its elements can be partitioned into the limit s and a remainder r n according to 

s n = s + r n , n e N . (2.8) 

A sequence transformation T is a rule which transforms a sequence {s n }^L into some other 
sequence {s' n }%L Q : 

r:{ Sn }- ^K}- . (2.9) 

A sequence transformation is only useful for our purposes if the transformed sequence {s' n }^L Q 
converges to the same limit s = Soo as the original sequence {s n }^L . If this is true, then the 
elements of the transformed sequence can also be partitioned into the limit s and a transformed 
remainder r' n according to 

s' n = s + r' n , n e No . (2.10) 

With the exception of some more or less trivial model problems, the transformed remainders 
r' n will be different from zero for all finite values of the index n. Thus, it looks as if nothing 
substantial is gained by applying a sequence transformation. However, the new remainders can 
have much better numerical properties than the original remainders. 

Accordingly, a sequence transformation is said to accelerate convergence if the transformed 
remainders {^j^Lo vanish more rapidly than the original remainders {^n}^Lo according to 

lim ^ = lim = 0. (2.11) 

n^oo r n n^oc s n — S 

In the case of some sequence transformations, it is actually possible to compute the transformed 
remainders {r' n }^— recursively if the numerical values of the original remainders {r n }^L are 
known Q. 

Thus, a sequence transformation tries to construct an approximation to the actual remainder 
r n and to eliminate it from s n . This yields a new sequence {s^jJJLo with hopefully superior nu- 
merical properties. In this context, we are confronted with the practical problem of constructing 
approximations to the actual remainders if we do not know explicit analytical expressions for 
them. However, even if we know explicit expressions, this normally does not help. For example, 
if the input data are the partial sums 

n 

Sn = Qk ' ( 2 ' 12 ) 

k=0 
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of an infinite scries, then the remainders are given by 

oo 

r n = - a k • ( 2 - 13 ) 

k=n+l 

We automatically have an explicit expression for the remainders if we have an explicit analytical 
expression for the terms a&. However, this does not help. The computation of r n via the infinite 
series ( |2.13| ) is normally not easier than the direct computation of the infinite series by adding up 
the terms successively. Consequently, the straightforward computation of the remainders of an 
infinite series is normally either not possible or not feasible, and we must use indirect approaches 
to obtain suitable approximations to the remainders. 

If we want to use a sequence transformation to speed up convergence, the minimal informa- 
tion, which we need to have, is a finite string of numerical values of input data, for example the 
sequence elements s n , s n +i, . . . , s n +fc with k,n £ No- This alone is, however, not enough for the 
construction of an approximation to the unknown truncation error, in particular since there can 
be an in principle unlimited variety of different types of remainders. Consequently, it should be 
obvious that a universal algorithm for the determination and elimination of completely arbitrary 
unknown remainders {r„}^L cannot exist. Delahaye and Germain-Bonne j7^,^o| were able to 
proof this rigorously. 

Fortunately, in most cases at least some structural and asymptotic information about the 
dependence of the remainders r n on the index n is available, or - if this is not the case - some 
more or less plausible assumptions can be made. We will show in the next sections how powerful 
sequence transformations can be constructed on the basis of simple assumptions. For example, 
let us assume that the elements of a sequence {s„}^ can be expressed by the following series 
expansion: 

oo 

s n = s + c j , n G N . (2.14) 

3=0 

Here, the coefficients Cj are assumed to be unknown, whereas the functions fj(n) are assumed 
to be known, but otherwise they are essentially arbitrary. 

A numerical algorithm can only involve a finite number of arithmetic operations. Conse- 
quently, a complete numerical determination of the remainder r n = ^2j^ Cj(pj(n) is impossible 
since this would require the determination of an infinite number of unknown coefficients Cj. So, 
the best we can hope for is the determination of a finite number of unknown coefficients Cj. 
If, however, the functions <Pj(n) decrease with increasing index j, we achieve an acceleration of 



convergence by eliminating the first k terms in the infinite series on the right-hand side of (2.14) 
which yields the transformed sequence 



oo 



s 'n = s + 22 4+ K fk+ K {n) , n e N . (2.15) 



K=0 



For example, if the functions <Pj(n) do not only decrease with increasing index, but are an 
asymptotic sequence according to 

<Pj+i(ri) = o(ipj(n)) , nGo, (2.16) 



then s n is according to ( 2.14 ) of order 0((/?o("-)) as n — > oo, whereas s' n is according to ( 2.15 ) 
of order 0(<^fe(n)), which for sufficiently large values of k leads to a substantial improvement of 
convergence. 
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The elimination of the leading k terms on the right-hand side of ( 2.14 ) can be accomplished 



with the help of a sequence transformation which is exact for the model sequence 

fc-i 

s n = s + Cjipj(n). (2-17) 

The elements of this model sequence contain k + 1 unknowns, the limit s and the k coefficients 
Co, C\, .. . , Cfe_i. Since all unknowns occur linearly, it follows from Cramer's rule that a sequence 
transformation can be constructed which is given as the ratio of two determinants, and which 



can determine the limit s in (2.17) if the numerical values of k + 1 sequence elements s„, s n +i, 
. . . , Sn+fc are available. An admittedly complicated recursive scheme for this general sequence 
transformation, which is determined completely as soon as the functions fj(n) are specified, was 
derived independently by Schneider Brezinski J82|, and Havie J83|]. 

Thus, model sequences s n — s + r n with remainders f n containing a finite number of terms as 



in ( 2.17 ) are very useful for the construction of sequence transformations that are special cases of 
the general sequence transformation mentioned above. If the remainders {f n }^Lo of the model 
sequence {s n }^o provide sufficiently accurate approximations to the remainders {r n }^ = Q of the 
sequence {s n }^Lo t° be transformed, then it can be hoped that a sequence transformation, which 
is exact for the model sequence {s n }^L , will be able to accelerate the convergence of {s n }^L 
to its limit s. More detailed discussions of this construction principle as well as many examples 
can be found in J6^,|7^|. 

Before we explicitly construct sequence transformations with the help of model sequences of 
the type of (|2.17 ), let us first mention that in the literature on sequence transformations, simple 



asymptotic conditions are used to classify the type of convergence of a sequence. For example, 
many practically relevant sequences {s n }^Lo converging to some limit s can be characterized by 
the asymptotic condition 

I im s " +1 ~ s = p , (2.18) 

n^oo s n — S 

which closely resembles the ratio test in the theory of infinite series. If \p\ < 1, the sequence is 
called linearly convergent, and if p = 1, it is called logarithmically convergent. 

3 Aitken Extrapolation and Wynn's Epsilon Algorithm 

Let us assume that the remainders of the elements of a sequence {s n }^L consist of a single 
exponential term: 

s n =s + cX n , c^O, |A|^1, neN . (3.1) 

For n — > co, this sequence obviously converges to its limit s if < |A| < 1, and it diverges away 
from its generalized limit s if |A| > 1. Thus, if this sequence converges, it converges linearly 
according to ( [2.18 ). Moreover, this model sequence contains the partial sums X)"=o(~ Z ) U = 



[1 — (— z)" +1 ]/[l + z] of the geometric series as special cases. 



If we consider s, c, and A in (3.1) as unknowns of the linear system s n+ j — s + cX J with 



j = 0, 1, 2, then we can easily construct a sequence transformation, which is able to determine the 



(generalized) limit s of the model sequence (3.1) from the numerical values of three consecutive 



sequence elements s n , s„+i and s n +2- A short calculation shows that 



A [ n) = s n - fiSNo, (3.2) 
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is able to determine the limit s of the model sequence fl3.l|) according to A^ = s. The forward 



difference operator A in ( |3.2[ ) is defined by its action on a function g = g(n): 

Ag(n) = g(n+l) - g(n) . (3.3) 



The A 2 formula ( |3.2| ) is certainly one of the oldest sequence transformations. It is usually 
attributed to Aitken ]5S{l , but it is actually much older. As already mentioned in Section [j], 
it was used in 1674 by Seki Kowa, the probably most famous Japanese mathematician of that 
period |56, pp. 90 - 91], and according to Todd |^4|, p. 5] it was in principle already known to 
Kummer | 8q| . 

The power and practical usefulness of Aitken's A formula is of course limited since it is 
designed to eliminate only a single exponential term from the elements of the model sequence 
(3.1). However, the quantities A x can again be used as input data in ( |3.2| ). Hence, the A 2 



process can be iterated, yielding the following nonlinear recursive scheme W% Eq. (5.1-15)] 



(«) 



A 



A {n) 



A 



(n) 



n £ N , 
A 2 4 n) 

i(") 



k, n £ N . 



(3.4a) 
(3.4b) 



In the case of doubly indexed quantities like A k , it will always be assumed that the difference 
operator A only acts on the superscript n but not on the subscript k: 



AA 



(n) 



A 



(n+l) 



A 



(n) 



(3.5) 



The numerical performance of Aitken's iterated A 2 process was studied in 74 , 36| . Concerning the 
theoretical properties of Aitken's iterated A 2 process, very little seems to be known. Hillion p7| 
was able to find a model sequence for which the iterated A 2 process is exact. He also derived a 
determinantal representation for -4.1 ■ However, Hillion's expressions for A^ contain explicitly 



the lower order transforms A^ , . . . , A^J_ X , . . . , A^^ , . . . , A^~ x k ' . Consequently, it seems that 



Hillion's result 



although interesting from a formal point of view - cannot help much to 



understand the theoretical properties of A^ . 

A more detailed discussion of Aitken's iterated A 2 process as well as additional references 
can for instance be found in Section 5 of (74|] or in j78| . The iteration of other sequence trans- 
formations is discussed in pal. 



An obvious generalization of the model sequence (3d) would be the following model sequence 
which contains k exponential terms: 



k-l 



s n — S + ^ Cj A T j 



n £ No . 



(3.6) 



Here, we assume that |Ao| > |Ai| > ... > |Afc_i| holds. Although the A 2 process (3^) is by 
construction exact for the model sequence (3.1), the iterated A 2 process (3.4) is not exact for the 
model sequence (3.6). Instead, this is true for Wynn's epsilon algorithm J6l|, which is another 
generalization of the A 2 process (3.2) and which corresponds to the following nonlinear recursive 
scheme: 



>) 
» 

z k+l 



o, 

(n+l) 
e fc-l 



» 



(n+l) 
k 



neNo 
»l 



k, n £ N . 



(3.7a) 
(3.7b) 
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(n) 

Only the elements e^ k with even subscripts provide approximations to the limit s of the sequence 
{s„}^Lq to be transformed. The elements e^fc+i with odd subscripts are only auxiliary quantities 
which diverge if the whole process converges. 

A straightforward calculation shows that A± — e\ ■ Accordingly, Aitken's iterated A 2 
process may also be viewed as an iteration of ■ However, for k > 1, A. k and e^u are m 
general different, although they of course have similar properties in convergence acceleration and 
summation processes. 

If the input data s n are the partial sums f n (z) = J2"=o ^ vzV °f ^ ne (formal) power series for 
some function f(z), s n = f n (z), then Wynn J6l| could show that his epsilon algorithm produces 
Pade approximants according to 



e 



(n) 
2k 



[n+k/k] f (z). (3.8) 



For Pade approximants, we use the notation of Baker and Graves-Morris p2^ , i.e., the Pade 
approximant [l/m]f{z) to some function f(z) is the ratio of two polynomials Pi{z) = pa + p\z + 
. . . + piz 1 and Q m (z) = 1 + q\Z + . . . + q m z m of degrees I and m in z according to 

[l/m] f (z) = Pi(z)/Q m (z). (3.9) 

Since the epsilon algorithm can be used for the computation of Pade approximants according 



to (3.8), it is discussed in books on Pade approximants such as |62j] . However, there is also 
an extensive literature which deals directly with the epsilon algorithm. On p. 120 of Wimps 
book p7[ | it is mentioned that over 50 articles on the epsilon algorithm were published by Wynn 
alone, and at least 30 articles by Brezinski. As a fairly complete source of references Wimp 
recommends Brezinski's first book ]63| . However, this book was published in 1977, and since 
then many more articles on the epsilon algorithm have been published. Moreover, the epsilon 
algorithm is not restricted to scalar sequences but can be generalized to cover for example vector 
sequences. A very recent review of these developments can be found in |89||. 



Aitken's iterated A 2 process (|3.4|) as well as Wynn's epsilon algorithm (3.7) are powerful 



accelerators for sequences which according to ( 2.18 ) converge linearly. However, contrary to a 



widespread misconception, the epsilon algorithm is not necessarily the most efficient convergence 
accelerator for linear convergence. In 1973 Levin |)0| introduced a sequence transformation 
which uses as input data not only a finite substring s„, s n +i, . . . , s n +fc of the sequence to be 
transformed, but also explicit estimates uj ni w n +i, . . . , uj n +k for the corresponding remainders 
r„, r„+i, . . . , r n +fc. The explicit incorporation of the information contained in the remainder 
estimates into the convergence acceleration or summation process makes Levin's sequence trans- 
formation j90j as well as closely related generalizations j74],|9l|] potentially very powerful. Some 
recent applications of Levin-type transformations plus additional references can for example be 
found in f|||9f|j|]. 

However, the use of explicit remainder estimates is not only the major advantage but also 
the major weakness of Levin-type transformations. If it is possible to find remainder estimates 
that are good approximations to the actual remainders of the sequence to be transformed, then 
experience indicates that Levin- type transformations can be extremely powerful. In contrast, 
Levin-type transformations will perform quite poorly if good remainder estimates cannot be 
found. 

If we want to extrapolate oligomer calculations, then we have to struggle not only with 
potentially large numerical inaccuracies in the input data but also with the fact that oligomers 
with a small number of repeat units may behave very irregularly. As a result large prediction 
errors can occur leading to poor performance of Levin-type transformations. 
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Wynn's epsilon algorithm (3/7) only uses the elements of the sequence to be transformed as 
input data and it is not possible to feed in any additional information into the transformation 
process. Consequently, the epsilon algorithm is not as powerful as Levin-type transformation 
can be under optimal conditions. On the other hand, there is no prediction error. Moreover, as 
for example discussed in fj\ Section 15.2], the epsilon algorithm is also remarkably insensitive 
to rounding errors and it can tolerate input data which either have a low relative accuracy or 
which behave in a comparatively irregular way. Thus, the real strength of the epsilon algorithm 
is its robustness. Due to its remarkable robustness, the epsilon algorithm is often able to produce 
meaningful and reliable results in situations in which other sequence transformations, which are 
in principle more powerful, fail. Experience indicates that Aitken's iterated A 2 process ([O]) is 
also less robust and more susceptible to rounding errors than Wynn's epsilon algorithm. 



4 Richardson Extrapolation, Wynn's Rho Algorithm, and 
Related matters 

In Section |3[ we discussed sequence transformations which work in the case of sequences that 



converge linearly according to (2.18). Typical examples of linearly convergent sequences are the 
partial sums of power series with a nonzero but finite radius of convergence. The computational 
problems, which occur in that case, can be studied via the Gaussian hypergeometric series 

3 F 1 {a,b;c i z) = £ z m . (4.1) 

„t- (c) m m! 

This series terminates if either a or 6 is a negative integer. Otherwise it converges linearly in the 
interior of the unit circle, i.e., for \z\ < 1, and it diverges for \z\ > 1. Moreover, the convergence 
of this series can for most values of the parameters a, 6, c G C be accelerated quite efficiently by 
sequence transformations [ p7| . 

Generally speaking, the acceleration of linear convergence is comparatively simple, both the- 



oretically and practically, as long as p in (2.18) is not too close to 1 (this would be a borderline 
case to logarithmic convergence). With the help of Germain-Bonne's formal theory of conver- 
gence acceleration |)8| and its extension [f74] , Section 12], it can be decided rigorously whether 
a sequence transformation is capable of accelerating linear convergence or not. Moreover, many 
sequence transformations are capable of accelerating linear convergence effectively. 

The situation is much more difficult in the case of logarithmic convergence. It follows from 



the expansion (2J3) for the truncation error that the partial sums X)"=o( I/ ~'~ -0~ z °f ^ ne Dirichlet 
series (2d) converge logarithmically to C,{z). As already discussed in Section |2[ the convergence 
of the Dirichlet series can become so slow that the evaluation of ((z) by successively adding up 
the terms (y + l)~ z is practically impossible. Of course, analogous problems can also occur in 
the case of other logarithmically convergent sequences and series, 

Moreover, there are also some principal theoretical problems. Delahaye and Germain-Bonne 
[^9[Q showed that no sequence transformation can exist which is able to accelerate the con- 
vergence of all logarithmically convergent sequences. Consequently, in the case of logarithmic 
convergence the success of a convergence acceleration process cannot be guaranteed unless ad- 
ditional information is available. Also, an analogue of Germain-Bonne's beautiful formal theory 
of the acceleration of linear convergence and its extension JTJ, Section 12] cannot exist. 

In spite of these complications, many sequence transformations are known which work rea- 
sonably well at least for suitably restricted subsets of the class of logarithmically convergent 
sequences. Examples, which will be discussed later, are Richardson extrapolation ]57|, Wynn's 
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rho algorithm |9!| and its iteration (74[ Section 6], as well as Osada's modification of the rho 
algorithm 100 1 , and the modification of the A 2 process by Bj0rstad, Dahlquist, and Grossc [101]. 



Nevertheless, there is a considerable amount of theoretical and empirical evidence that sequence 
transformations speed up logarithmic convergence in general less efficiently than linear conver- 
gence. 

Another problem with logarithmic convergence is that numerical instabilities are more likely 
than in the case of linear convergence. A sequence transformation can only accelerate convergence 
if it succeeds in extracting some information about the index-dependence of the truncation errors 
from a finite set s n , s n +i, . . . , s n +k of input data. Normally, this is done by forming arithmetic 
expressions involving higher weighted differences. However, forming higher weighted differences 
is a potentially unstable process which can easily lead to a serious loss of significant digits or even 
to completely nonsensical results. If the input data are the partial sums of a strictly alternating 
series, the formation of higher weighted differences is normally a remarkably stable process. If, 
however, the input data all have the same sign, numerical instabilities due to cancellation are 
quite likely, in particular if convergence is very slow. Thus, if the sequence to be transformed 
converges logarithmically, numerical instabilities are a serious problem and at least some loss of 
significant digits is to be expected. 

In p5| , a combined nonlinear-condensation transformation was described which first converts 
a slowly convergent monotone series into a strictly alternating series. Since the subsequent 
transformation of a strictly alternating series is a remarkable stable process, it is in this way 
possible to evaluate special functions, that are defined by extremely slowly convergent series, not 
only relatively efficiently but also very close to machine accuracy. Unfortunately, an analogous 
approach is not possible if oligomer calculations produce logarithmically convergent sequences: 
The conversion of the monotone series to the alternating series requires that series terms with 
large indices can be computed. This is of course not possible in the case of oligomer calculations. 

For the construction of sequence transformations, the standard interpolation and extrapola- 
tion methods of numerical mathematics are quite helpful. Let us assume that the values of a 
function f(x) are only known at some discrete points xo < x\ < ■ ■ ■ < x m , and that we want 
to estimate the value of / at some point £ ^ {xo, X\, . . . , x m }. If xq < £ < x m , this problem is 
called interpolation, and if either £ < xq or x m < £, this problem is called extrapolation. These 
topics are discussed in any book on numerical analysis. More specialized treatments can for 



instance be found in [102, 103 1 



For the construction of extrapolation processes, we postulate the existence of a function S of 
a continuous variable which coincides on a set of discrete arguments {x n }^L with the elements 
of the sequence {s n }^_ to be transformed: 

S{x n ) = s n , neN . (4.2) 

This ansatz reduces the problem of accelerating the convergence of a sequence to an extrapolation 
problem. If a finite string s m ,s m +i, ■ ■ . , s m +fc of fc + 1 sequence elements is known one can 
construct an approximation Sk(x) to <S(x) which satisfies the k + 1 interpolation conditions 

S k (x m+j ) = s m+ j, 0<j<k. (4.3) 

In the next step, the value of Sk(x) has to be determined for x — > Xoq. If this can be done and 
if there exists a function S(x) which can be approximated at least locally by a suitable set of 
interpolating functions, one can expect that Skfeoo) will provide a better approximation to the 
limit s — of the sequence {s„}^ than the last sequence element s m +fe which was used for 
the construction of Sk(x). 
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The most common interpolating functions are either polynomials or rational functions. These 
two sets will also lead to different convergence acceleration methods. If a sequence transformation 
is constructed on the basis of polynomial interpolation, it is implicitly assumed that the fc-th 
order approximant Sk(x) is a polynomial of degree k in x: 

S k {x) = c + ax + ■ ■ ■ + c k x k . (4.4) 

For polynomials, the most natural extrapolation point is x — 0. Consequently, the interpolation 
points x n have to satisfy the conditions 

x Q > xi > ■ ■ ■ > x m > x m+ i > • ■ ■ > , (4.5a) 
lim x n = 0. (4.5b) 



n — 'oc 



The choice Xoo = as the extrapolation point implies that the approximation to the limit is to 



be identified with the constant term Co of the polynomial (4.4). 

Several different methods for the construction of interpolating polynomials Sk(x) are known 
in the mathematical literature. Since only the constant term of a polynomial Sk(x) has to be 
computed and since in most applications it is desirable to compute simultaneously a whole string 



of approximants 50(0)551(0), £2(0), . . . , the most economical choice is Neville's scheme [ 104 1 for 
the recursive computation of interpolating polynomials. If wc set x — 0, Neville's algorithm 
reduces to the following 2-dimensional linear recursive scheme ^J, p. 6]: 

./V (n) (s n ,x n ) = s n , neN a , (4.6a) 
.,(„), v x n Af( n+1 \s n+1 ,x n+1 ) - x n+k +iJ^k (s n ,x n ) 

N^iSntXn) = - , fc,neN . (4.6b) 

Xn Xn+k+1 

In the literature on convergence acceleration this variant of Neville's recursive scheme is usually 
called Richardson extrapolation J57J. Obviously, A/^ (s„, x n ) is exact if elements s n of the 
sequence to be transformed are polynomials of degree k in the interpolation points x n : 

k-l 

S n = s + Y1 C 3 X n +1 ' fc '" G N ° ' ( 47 ) 

3=0 

As is well known, some functions can be approximated by polynomials only quite poorly, but 
by rational functions they can be approximated very well. Consequently, it is to be expected 
that at least for some sequences {s n }^Lo ra tional extrapolation will give better results than 
polynomial extrapolation. Let us therefore assume that the approximant Sk(x) can be written 
as the ratio of two polynomials of degrees I and to, respectively: 

a + a\x + a 2 x 2 H h aix 1 

Sk{x) = b + b lX + b 2X > + ... + b mX r» > M,m6N . (4.8) 

This rational function contains I + m + 2 coefficients ao, . . . , a; and 60, . . . , b m . However, only 
l+m+1 coefficients are independent since they are determined only up to a common nonvanishing 
factor. Usually, one requires either 60 = 1 or b rn = 1. Consequently, the k + 1 interpolation 



conditions (4.3) will determine the coefficients ao, . . . , a; and 60, . . . , b m provided that k = I + m 
holds. 

The extrapolation point x^ = is also the most obvious choice in the case of rational 
extrapolation. Extrapolation to x^ = implies that the interpolation points {£ ra }^ m (4 



have to satisfy (Eq) and that the approximation to the limit is to be identified with the ratio 



ao/bo of the constant terms of the polynomials in (4.8) 
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If I = m holds in (4.8), extrapolation to infinity is also possible. In that case the interpolation 



points {2 n }^Lo have to satisfy 



< x < xi < ■ ■ ■ < x m < x m+ i < • • • , (4.9a) 
lim x n = oo . (4.9b) 



n — >oo 



In the case of extrapolation to infinity, only the coefficients a; and bi of the polynomials in (4.8), 
that are proportional to highest power x , contribute. Consequently, the approximation to the 
limit has to be identified with the ratio ai/bi. 

As in the case of polynomial interpolation, several different algorithms for the computation 
of rational interpolants are known. A discussion of the relative merits of these algorithms as well 



as a survey of the relevant literature can be found in [105, Chapter III]. 

The most frequently used rational extrapolation technique is probably Wynn's rho algorithm 
©: 

p^l = 0, = s n , neN , (4.10a) 

(n) , Xn+k+1 — X n . . , , 

Pk+i = Pk-i + "(n+ii M> k,neN . (4.10b) 

Pk ~ Pi 



Formally, the only difference between Wynn's epsilon algorithm (3.7) and Wynn's rho algorithm 
is that the rho algorithm involves a sequence of interpolation points {xn}^^ which have to 
satisfy (4.£). As in the case of the epsilon algorithm, only the elements p^ with even subscripts 



serve as approximations to the limit. The elements p^fc+i with odd subscripts are only auxiliary 
quantities which diverge if the whole process converges. 

Wynn's rho algorithm is designed to compute even-order convergents of Thiele's interpolating 
continued fraction |l06|] and to extrapolate them to infinity. The even-order convergents are 
rational functions of the following type: 

c i\ afcX fc + ai~-ix k ~ 1 H h a x x + a Q 

bkX K + b k ^ix K 1 + • • • + bix + b a 

Thus, the ratio a^jb^ is to be identified with the approximation to the limit. 

The epsilon and the rho algorithm have complementary features. The epsilon algorithm is a 
powerful accelerator for linear convergence and is also able to sum many divergent alternating 
series, whereas the rho algorithm fails to accelerate linear convergence and is not able to sum 
divergent series. However, it is a very powerful accelerator for some logarithmically convergent 
sequences. 

The properties of Wynn's rho algorithm are for example discussed in books by Brezinski 
(see |H pp. 102 - 106] and [f| pp. 96 - 102] and Wimp ||, pp. 168 - 169]. In these books the 
connection of the p algorithm with interpolating continued fractions is emphasized and it is also 

(n) 

shown that the transforms p 2t J can be represented as the ratio of two determinants. Moreover, 



there is an article by Osada [107], discussing its convergence properties. But otherwise, relatively 



little seems to be known about this sequence transformation. 

In the case of the rho algorithm, we can proceed as in the case of the Aitken formula (3.2) 
and construct an iterated transformation. For k = 1, we obtain from ( 4.10] ) : 



(x n +2 ~ x n )[As n+1 ][As n ] 

P\ = S n+1 + TT TTT 1 FT TTT T , Tl G N . 4.12) 

[Ax n+1 \[As n \ - [Ax n \[As n+1 \ 



This expression, which can be viewed to be a kind of weighted A 2 formula, can be iterated 
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yielding |88|, Section 6.3] 



(n) 



n G N . 



(4.13a) 



(n) 
fc+1 



(n+1) 



fc, 7i G N . 



(x n+2k+2 - x n )[AWi n+1) ][Awi n) ] 



{x n +2k+2 - X n+ i)[AW i k n ' > ] - {x n +2k+l - X n ) [AW, 



(n+1) 
k 



» 



However, in 



(4.13b) 

it was shown that other 



This is not the only possibility of iterating p 2 
iterations of - for instance the one derived by Bhowmick, Bhattacharya, and Roy [108] 
are significantly less efficient than Wl , which has similar properties as Wynn's rho algorithm 
(4.10), from which it was derived. 

The main problem with the practical application of those sequence transformations, that are 
based upon interpolation theory, is that one has to find a sequence {^nl^o of interpolation 
points which produces good extrapolation results for a given sequence {s n }^Lo of input data. 
For example, in most practical applications the Richardson extrapolation scheme (^1]) is used in 
combination with the interpolation points x n = 1/ (n + 0), where (3 is a positive shift parameter. 
Then, AQ (s n , l/(n + (3)) possesses a closed form expression (see for example 68, Lemma 2.1 
on p. 313] or @, Eq. (7.3-20)]) 



Aj>*>(/?,* n ) = Mi n \s n ,l/{n + (3)) 



3=0 



>n+j 



and the recursive scheme (4.6) assumes the following form 

A n) (/3,s„) = s„, n G No , 

A^ (/?,*„) 



fc,nGN , 
Eq. (7.3-21)] 



(4.14) 



(4.15a) 

A(" +1) (/3,.s„ +1 ) + ^±^AA( n) (/3, Sn ), fc,nGN . (4.15b) 

Similarly, Wynn's rho algorithm ( 4. 10] ) and its iteration (4.13) are normally used in combination 
with the interpolation points x n = n + 1 , yielding the standard forms (see for example ] [74| Eq. 
(6.2-4)]) 



(n) 
P-l 

(n) 
Pk+1 



0, 

(n+1) 
Pk-1 



(n) 

Po 



k + l 



n G N , 



(n+1) 
Pk 



(«) 



k, n G N . 



and J88|, Section 6.3] 



(n) 



n G N . 



(n) 
fe+1 



(n+1) 



(2k + 2)[AWl n+i> ][AW, 



;(»)! 



fc, n G N . 



(4.16a) 
(4.16b) 

(4.17a) 
(4.17b) 



(2k + 1)A 2 W^ 

In the literature on convergence acceleration and extrapolation, Richardson extrapolation 
as well as Wynn's rho algorithm and its iteration are normally used in their standard forms 
( 4.15 ), ( 1.16 ), and ( 4.17 ). However, it is important to recognize that their general forms (4.6), 
( 4.10 ), and ( 4.13 ), respectively, which all involve unspecified interpolation points, may be required 
for oligomer calculations. Let us assume that the elements s n of a logarithmically convergent 
sequence can be expressed by series expansions of the following kind, 



s + (n + (3)- a J2^/(n + PY 



n G N . 



(4.18) 



Extrapolation Methods for Oligomer Calculations 



17 



where a is a positive decay parameter, [3 is a positive shift parameter, and the Cj are unspecified 
coefficients. In general, we have no a priori reason to assume that oligomer calculations would 
only produce sequences of that kind with an integral decay parameter, and we want to accelerate 
the convergence of sequences of that kind even if a is nonintegral. 



In Theorem 14-4 of |74|, it was shown rigorously that the standard form ( 4.15 ) of Richardson 
extrapolation is able to accelerate the convergence of the model sequence ( 4.18 ) if a = 1, 2, ... , 



but fails if a is nonintegral. In the case of Wynn's rho algorithm and its iteration, we are not 
aware of any explicit proof, but there is considerable experimental evidence that the standard 



forms ( 4.16 ) and ( 4.17 ) of these transformations also fail to accelerate the convergence of the 
model sequence ( |1.18| ) if the decay parameter a is nonintegral (see Section 14.4 of |7^| ). 



If the decay parameter a of a sequence of the type of (4.18) is known, then Osada's variant 



of Wynn's rho algorithm can be used | IOC , Eq. (3.1)] 



p^l = 0. 



-(«) 

Pk+i = 



Pt\ 1] 



_(n) _ 

k + a 



n G N , 



-(n+l) 
Pk 



-(«) 
Pk 



k, n 6 N . 



(4.19a) 
(4.19b) 



For a = 1, Osada's variant becomes identical with the standard form (1.16). Osada also demon- 



strated that his variant accelerates the convergence of the model sequence (4.18), and that the 

(4.20) 



transformation error satisfies the following asymptotic estimate [100, Theorem 4]: 

= 0(n- a ~ 2k ) , n -> oo. 



_(n) 

p\k 



Osada's variant of the rho algorithm can also be iterated. From ( 4.19 ) we obtain the following 



expression for in terms of s„, s n +i, and s n+2 : 



-in) 



>n+l 



(a+1) [As n ][As n+1 ] 



[A 2 s„] 



n G N . 



(4.21) 



If this iteration is done in such a way that a is increased by 2 with every recursive step, we 
obtain the following recursive scheme [pq, Eq. (2.29)] which was derived originally by Bj0rstad, 



Dahlquist, and Grosse [101, Eq. (2.4)] 



(n) 




n G N . 



W (n) _ W in+1) (2fc + a + l) [Awi' l+1) ][Awi n) ] 
vv k+l — VV k 



(2k + a) 



A 2 Wi n) 



k, n G N . 



(4.22a) 
(4.22b) 



For a = 1, this weighted iterated A 2 algorithm becomes identical with the standard form (4.17) 
of the iteration of Wynn's rho algorithm. 

( n ) 

Bj0rstad, Dahlquist, and Grosse showed that W k accelerates the convergence of the model 



sequence ( [4.18[ ), and that the transformation error satisfies the following asymptotic estimate 

(4.23) 



101, Eq. (3.1)] 



w[ n) 



= 0(n- a - 2k ) 



n — > oo . 



The two sequence transformations and W^' w require the same set s n , s„+i, ... , s n +2fe of 
elements of the model sequence ( 4.18| ) as input data. Since the asymptotic error estimates (4.20) 
and (4.23) are equivalent, the two transformations are also asymptotically equivalent, i.e., in the 
limit of large indices n, although they of course may produce different results for small indices 
n. Moreover, these two transformations are also optimal, since the elimination of two terms of 



the infinite series on the right-hand side of (4.15) requires the input of the numerical values of 
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two sequence elements. Thus, no sequence transformation, which only uses the numerical values 
of the elements of the model sequence ( 1.18| ) as input data, could improve on the asymptotic 
(n — > oo) error estimates ( 4.2CQ and ( 4.23| ). 

The knowledge of the decay parameter a is crucial for an application of the transformations 



( 4.19 ) and ( 4.22 ) to a sequence of the type of ( 4. IS ). An approximation to a can be obtained 
with the help of the following nonlinear transformation: 



T„ = 



[A 2 s n ] [A 2 s r 



[A« n+ i] [A 2 s„ +1 ] - [As„ +2 ] [A 2 s„] 



1, 



n € N 



(4.24) 



The transformation T n , which was first derived in a somewhat disguised form by Drummond [|109j 
and later rederived by Bj0rstad, Dahlquist, and Grosse [101|, is essentially a weighted A 3 method. 



Consequently, it is potentially a very unstable method if the relative accuracy of the input data 



is low, and stability problems can never be excluded. Bj0rstad, Dahlquist, and Grosse [101, Eq. 
(4.1)] also showed that 



a = T n + 0{l/n 2 ) 



(4.25) 



if the elements of the model sequence (4.18) are used as input data 



5 An Example: Extrapolation of Oligomer Calculations 
for Unit Cell Energy of Polyacetylene 

All-trans polyacetylene with the repeat unit — (CH=CH) — , which contains 14 electrons, is the 
simplest gMasz-onedimensional polymer containing conjugated double bonds. As a result of the 
conjugation and the pronounced anisotropy, polyacetylene exhibits a variety of interesting proper- 
ties. Included amongst these properties is an enormous increase in conductivity induced by dop- 
ing (see for example Jl|, |ll0| ) and a very large electric field induced polarization [^],f47| compared 
to the corresponding non-conjugated polymer polyethylene with the repeat unit — (CH2-CH2)-. 
Accordingly, there is an enormous number of articles dealing with polyacetylene. A reasonably 
complete bibliography would clearly be beyond the scope of this article. 

In this article, we want to demonstrate that sequence transformations are useful numerical 
tools for the extrapolation of oligomer calculations. Consequently, we are not so much interested 
in obtaining highly accurate results. Instead, we want to study the typical behavior of quantum 
chemical calculations on a sequence of oligomers with increasing size. For this purpose, it is 
sufficient to use relatively small basis sets which allowed us to carry out calculations on fairly 
large oligomers expeditiously. 

Here we some consider calculations for the ground state energy of the oligomers H— (CH=CH)n— 
H. These will illustrate how new and improved results can be obtained but, by no means exhaust 
the range of possibilities. Other gMasi-onedimensional polymers, as well as the application to 
other properties such as hyperpolarizabilities and electronic transition energies, will be discussed 



elsewhere [111|. The calculated total ground state energy En of an oligomer consisting of N 
repeat units is an extensive quantity. This means that we first have to convert En to an intensive 
quantity before we can do extrapolations to the infinite chain limit. The most obvious approach 
consists in defining average energies or energies per repeat unit according to 

E<fi v) = E N /N. (5.1) 



If we compare the total energy of an oligomer H— (CH=CH)n — H with the total energy of a 
hypothetical ring molecule consisting of N repeat units — (CH=CH) — , and if we also assume 
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that TV is so large that angular strain can be neglected, then we see that the saturation of the 
terminal free valences by hydrogen atoms should lead to an energy contribution which is nearly 
constant, i.e., which depends on TV only quite weakly. Consequently, we may reasonably assume 
that the average energy E^ v ^ will contain a contribution which is proportional to 1/N arising 
from the chain ends and, therefore, can be expressed as 



+ C/N + $ 



(av) 
N 



(5.2) 



where C is an TV-independent constant and $jy is an unknown function of N which should 
decay faster than 1/TV as TV becomes large. 

If we want to extrapolate the average energies successfully to the infinite chain limit, then we 
must have an idea of how the unknown function depends upon TV. It is a relatively obvious 



idea to express 



or equivalently E^' by a power series in 1 /TV according to 



E 



(av) 
JV 



= E {a 



E 



(5.3) 



If this is case, then the sequence of the average energies converges logarithmically to the infi- 
nite chain limit according to (2.18). As discussed in Section 0, the acceleration of logarithmic 



convergence is normally more difficult than the acceleration of linear convergence. Nevertheless, 



the standard forms ( 4.15 ), ( 4.16 ), and ( 4.17 ) of Richardson extrapolation, Wynn's rho algorithm 



and its iteration should be effective accelerators for sequences of the kind of (5.3). 

However, there are still some unresolved theoretical questions. According to Cioslowski and 



Lepetit 112, p. 3544], who performed a perturbation theoretic analysis of the TV-dependence of 



Hartree-Fock oligomer energies, the power series expansion (5.3) for the average energy can 



not only contain powers of 1/N but also nonanalytic terms of the type exp(— 7 TV) with 7 > 0. 
For sufficiently large values of TV, these nonanalytic contributions will be negligible. However, 
they could well have a negative effect on the performance of convergence acceleration processes. 
Sequence transformations designed for the acceleration of logarithmic convergence are usually 
not good at eliminating nonanalytic contributions typical of linear convergence, and vice versa. 
Alternatively, Cui, Kertesz, and Jiang ElJ proposed to compute the energy per repeat unit 



not via the average energies (5.1) but via the energy differences 



E 



(dif) 

N 



= E 



_x -E N = {N + 1)E%% -NE™ 



(5.4) 



The argument of Cui, Kertesz, and Jiang was that the energies of two oligomers containing iV+1 
or TV monomer units, respectively, should essentially differ by the energy of an inner monomer 
unit 11S]. Equivalently, this means that the effects due to the finite size of the oligomers as well 



as the 1/iV errors due to the saturation of free terminal valences should depend on TV only quite 
weakly and are largely canceled by forming differences. 
If we now combine (5.2) with (|5.4|), we obtain 



E 



(dif) 

N 



E {a 



A[TV$ 



(av)] 
N J 



(5.5) 



where 3?^ is the function defined in (5.2) that describes the higher order effects of the TV- 



dependence of the average energies. Comparing (t 



(av) 



N has been assumed to decay faster than 1/TV, we see that E N 

(av) 



and (5.?:) and taking into account that 
" does not contain a 1/TV 



contribution and converges more rapidly than E N to the infinite chain limit. If, for example 



it is assumed that $^ vj ~ 1/TV 2 as TV —> 00, we obtain A [TV$^ V) ] ~ 1/TV 2 which shows that the 
energy differences converge faster than the average energies. 
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In order to illustrate some of the above procedures we carried out Hartree-Fock ground state 
energy calculations for the oligomers H— (CH=CH)n — H with N < 16 using an ST0-3G basis 
set. A fixed geometry was used for all oligomers with Rc=c = 1-366 A, Rc-c = 1-450 A, 
i?c_ H = 1-085 A, Z c =c-c = 123.90 deg, Z H -c-h = 117.13 deg, and Z H -c-h = 117.13 deg. 
This geometry implies Zh-c=c — 118.97 deg, except for the terminal hydrogen on cither end, 
for which Zh-c=c = 123.90 deg. In Table I, we list the total energies En, the average energies 
i?^ v ' defined in (5.1), and the energy differences E^ defined in (5.4). The data displayed in 
Table I show that the energy differences converge indeed much more rapidly than the average 
energies. 



Table I: Total Hartree-Fock energies En, average energies E 
and energy differences E^ 1 ^ using an ST0-3G basis set. 
For geometry see the text. 



N 


Em 




p(av) 
J^N 




p(dif) 


1 


-77.0672438490 


-77 


0672438490 


-75 


943944441 


2 


-153.011188290 


-76 


5055941450 


-75 


945112888 


3 


-228.956301178 


-76 


3187670593 


-75 


945528271 


4 


-304.901829449 


-76 


2254573623 


-75 


945641947 


5 


-380.847471396 


-76 


1694942792 


-75 


945676982 


6 


-456.793148378 


-76 


1321913963 


-75 


945688518 


7 


-532.738836896 


-76 


1055481280 


-75 


945692475 


8 


-608.684529371 


-76 


0855661714 


-75 


945693869 


9 


-684.630223240 


-76 


0700248044 


-75 


945694368 


10 


-760.575917608 


-76 


0575917608 


-75 


945694549 


11 


-836.521612157 


-76 


0474192870 


-75 


945694615 


12 


-912.467306772 


-76 


0389422310 


-75 


945694639 


13 


-988.413001411 


-76 


0317693393 


-75 


945694649 


14 


-1064.35869606 


-76 


0256211471 


-75 


945694650 


15 


-1140.30439071 


-76 


0202927140 


-75 


945694650 


16 


-1216.25008536 


-76 


0156303350 







In the next step, we want to analyze whether the average energies and the energy differences 
behave like the elements of the model sequence ( 4.1S ), i.e., whether we can identify a well- 
defined decay parameter a. For that purpose, we use the transformation T n defined in (4.24). 
All sequence transformations discussed in this article are gauged in such a way that the indices 
of the elements of the sequence to be transformed start with zero. Thus, we use the average 
energies and the energy differences as input data for T„ according to either s n — E^7\ with 
< n < 16 or s n = E^ with < n < 14. 

The results in the third column of Table II show quite clearly that the average energies 
possess a contribution which is proportional to \/N. On the basis of these results alone, one 
might be tempted to generalize this observation and conclude that the average energies can be 
expressed by expansions in powers of l/N according to (5.3). However, the approximations to 
the decay parameter of the energy differences in the fifth column show an erratic behavior and 
there is no indication that they might converge. Thus, we have no reason to assume that energy 



Extrapolation Methods for Oligomer Calculations 



21 



Table II: 
energies 



Approximations to the decay parameter a of the average 
E„7\ and the energy differences E^+l according to ( 4.24 ) 



n 




p(av) 








p(dif) 

- C/ n+l 




u 


— It 


00 ( z4oo4y0 


i 


UUZ00Z4 


ik 
— to 


CiA QO/1 A A A 1 

y4oy44441 


— 0. ( zUool < 


1 


— 10 


rnr rO/1 i /l rn 

o0ooy414o0 


U 


0070070 

yy i zv i\) 


IK 

— to 


O/I^IIOOOO 


1 Q C/IOQIQ 
lo.04y<3l<3 


z 


— 10 


olo ( D ( Uoyo 


o 
U 


yy to luz 


IK 

— to 


O/i r;f:0Q07i 


oi noon7r; 
zl.UzzU(i3 


Q 

o 


— 10 


nop: /l r;7QSOQ 

ZZ040 ( oozo 


o 
U 


OOC A 1 0£ 

yyo4iut> 


— 1 


O/l KP.A 1(1/17 

y4oo4iy4 < 


Q1 (\CiKCi r iCi 

ol.Uooooo 


4 


7R 

— to 


1 £0/10/10*700 


o 
U 


onono a 1 

yyyuz4i 


7K 

— 10 


O/l r:fi7finQO 

y4ob my»z 


A A QQKKCIO 

44.000oyz 


5 


-76 


1321913963 





9994399 


-75 


945688518 


72.270674 


G 


-76 


1055481280 





9996933 


-75 


945692475 


84.907033 


7 


-76 


0855661714 





9998391 


-75 


945693869 


210.38728 


8 


-76 


0700248044 





9999177 


-75 


945694368 


-403.50000 


9 


-76 


0575917608 





9999589 


-75 


945694549 


6.0000000 


10 


-76 


0474192870 





9999827 


-75 


945694615 


-2.6578947 


11 


-76 


0389422310 





9999829 


-75 


945694639 


-10.000000 


12 


-76 


0317693393 





9999976 


-75 


945694649 




13 


-76 


0256211471 






-75 


945694650 




14 


-76 


0202927140 






-75 


945694650 




15 


-76 


0156303350 













differences decay like a fixed power of 1/N. In view of and (5J5) this observation implies 



that the unknown function defined in (5J2), which describes the higher order effects of 



the truncation errors in the average energies, also does not decay like a fixed power of 1/N. 
Consequently, our observation contradicts the perturbation theoretic analysis of Cioslowski and 



Lepetit [112] who concluded that the average energies can be expressed by a power series in 1/N 



according to (5.3) plus possibly some nonanalytic contributions. 

If we want to apply sequence transformations to speed up the convergence of oligomer calcu- 
lations to the infinite chain limit, we have to assume that the iV-dependence of the truncation 
errors in the energy differences possesses a structure which is sufficiently regular to be detected 
and utilized in convergence acceleration processes. The energy differences in the last column 
of Table I showed a remarkably rapid convergence, indicating that the truncation errors might 
converge exponentially. A rigorous theoretical proof of this assumption is of course not in sight. 
So, we have to provide numerical evidence for our conjecture. 

There is a simple ratio test related to the ratio test for infinite series, which can provide 
evidence for exponential decay. If we form the ratio 

Sn+2 — Sn+1 ^ S 7i+1 



As r< 



(5.6) 



of the elements of the model sequence (31), then we obtain lZ n = A. If the input data do not have 
the simple form of the elements of this model sequence, but nevertheless decay exponentially, 
then we only get a constant as n — > oo. In Table III, we apply the ratio test (5.6) to the 
input data s„ = with < n < 11. The results in the last column of Table III seem 

to converge, albeit slowly. Therefore, they provide evidence that the energy differences indeed 
converge exponentially to the infinite chain limit. 
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Table III: Ratio test (5.6) of the energy differences 



n 


p(dif) 
nj n+ 1 


n n 


n 
u 


7^ Q/l 3Q/I/I/I/1 1 


U.oOOO 


i 
i 


7K Q/1 ^1 1 9888 


97*37 


2 


-75.945528271 


0.3082 


3 


-75.945641947 


0.3293 


4 


-75.945676982 


0.3430 


5 


-75.945688518 


0.3523 


6 


-75.945692475 


0.3580 


7 


-75.945693869 


0.3627 


8 


-75.945694368 


0.3646 


9 


-75.945694549 


0.3636 


10 


-75.945694615 




11 


-75.945694639 





There is another approach which can provide evidence that the truncation errors of the energy 
differences indeed decay exponentially. If we apply sequence transformations, that are known to 
accelerate linear convergence but fail in the case of logarithmic convergence like Aitken's iterated 



A process (3^) or Wynn's epsilon algorithm (3/7), and observe an acceleration of convergence 



then we have very strong evidence for our conjecture. 

There is a problem with the acceleration of the convergence of the energy differences. The 
last column of Table I shows that the energy differences E^ 1 ^ have for N = 14 already converged 
to all digits given. Thus, there are no additional digits left which could provide information to 
the extrapolation methods. All sequence transformations form higher weighted differences of the 
input data, and so does Aitken's iterated A 2 process or Wynn's epsilon algorithm. If, however, 
the input data have almost reached convergence and there are no additional digits left, numerical 
instabilities or even completely nonsensical results can never be excluded. 

If the input data so, si, ... , S2fc are known, then the transform with the highest possible 
subscript, that can be computed with the help Wynn's epsilon algorithm ( p.7[ ), is e^. . In general, 
this is the most efficient approach and one can normally expect the best results. If, however, the 
input data behave like the energy differences and differ at most with respect to a few digits, it is 
recommended that one apply the epsilon algorithm differently Instead, we start from the initial 
string el 2 -* = E^l with < n < n max and compute the string eff with < n < n max — 2. Then 
we analyze the convergence and reliability of this string and compute, if necessary and possible, 
the next string en with < n < 7i max — 4. This process can be repeated until either convergence 
has been achieved or the accumulation of rounding errors makes it impossible to proceed. 



In Table IV, we show the effect of Wynn's epsilon algorithm (3/?) on the energy differences 



The extrapolation results provide strong evidence that the truncation errors of the energy dif- 
ferences converge exponentially. For example, eg = —75.945694655, which requires the initial 
values Cq = = E^ dli \ . . . , e[, = E^ dli ^ for its computation, reproduces all but the 

last digit of eg 13 "* = E^ 1 ^ = —75.945694650. Thus, the epsilon algorithm makes it unnecessary to 
compute the oligomer energies En with 9 < N < 14, which is indeed a remarkable achievement. 
We also used Aitken's iterated A 2 process (3.4) for the extrapolation of the energy differ- 
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Table IV: Epsilon extrapolation of the energy differences. 



n 


(n) _ p(dif) 
fc — 




>) 

fc 2 




>) 


>) 

fc 6 


n 

U 


-75 


943944441 


— 1 


y^o ( o ( oyz 


7^ 


QzL^RQl ^97 


7^ QzL^RQAR^l 

— t o.y^ouy^iooi 


i 
i 


-75 


945112888 


7^ 


QzL^fi8zL777 


7^ 


y^rOoy^oiz 


7^ QzL^fiQzLfi^ 

— i o.y^ooy^ooo 


9 


-75 


945528271 


7^ 


y^ooyzoyu 


7^ 


QAKP.QAp.'lA 

y^rOoy^rUo^: 


7^ QA^fiQAK^l 

— i o.y^ooy^ooi 


9 
O 


— (5 


94564194/ 


— 1 


y^ooy^ioi 


7^ 


QzL^RQzLfi^9 




4 


-75 


945676982 


-75 


945694541 


-75 


945694651 


-75.945694653 


5 


-75 


945688518 


-75 


945694627 


-75 


945694654 


-75.945694652 


G 


-75 


945692475 


-75 


945694646 


-75 


945694653 


-75.945694653 


7 


-75 


945693869 


-75 


945694652 


-75 


945694653 


-75.945694763 


8 


-75 


945694368 


-75 


945694653 


-75 


945694653 




9 


-75 


945694549 


-75 


945694653 


-75 


945694654 




10 


-75 


945694615 


-75 


945694656 








11 


-75 


945694639 


-75 


945694650 








12 


-75 


945694649 












13 


-75 


945694650 













ences. These extrapolation results were virtually identical with the epsilon extrapolation results 
presented in Table IV. 

When we performed oligomer calculations for polyacetylene using more extensive basis sets 
and other geometries, or for other gwasi-onedimensional polymers, we again observed an expo- 
nential decay of the energy differences [111]. Thus, this exponential decay seems to be a general 
feature which is not restricted to polyacetylene. 



6 Summary and Conclusions 

(Jttasi-onedimensional stereoregular polymers like polyacetylene are currently of considerable 
interest, not only scientifically but also because of numerous possible technological applications. 
There arc basically two different approaches for doing ab initio electronic structure calculations 
on such systems which have largely complementary advantages and disadvantages. One method, 
the so-called crystal orbital method, is based on the concepts of solid state theory. It uses periodic 
boundary conditions and leads to a band structure treatment. Accordingly, this approach is well 
suited to describe those features that depend crucially on the in principle infinite extension of 
the polymer. However, there are still some open computational problems with the crystal orbital 
method. For example, it is very difficult to properly take into account the effect of nonperiodic 
perturbations such as the interaction with an electric field or to incorporate high level electron 
correlation effects. 

The other method is essentially a quantum chemical approach. It approximates the polymer 
by oligomers consisting of a finite number of monomer units, i.e., by molecules of finite size. In 
this way, the highly developed technology of quantum chemical molecular ab initio programs can 
be used. Thus, electron correlation can be treated using any of the techniques available in these 
programs. In addition, nonlinear optical properties may be calculated via this approach. At the 
present time, this cannot be done reliably within the framework of crystal orbital theory even at 
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the Hartree-Fock level. 

Unfortunately, oligomers of finite size are not necessarily able to model those features of 
a polymer which crucially depend on its in principle infinite extension since finite oligomer 
calculations can suffer badly from slow convergence to the infinite chain limit. However, one can 
perform electronic structure calculations for a sequence of oligomers with an increasing number 
of monomer units and try to determine the limit of this sequence with the help of suitable 
extrapolation methods. 

Many different approaches for doing the extrapolations are possible. In this article, we use 
sequence transformations. The mechanism, whereby sequence transformations accomplish an 
acceleration of convergence, can be understood via the Euler-Maclaurin formula which is able to 
produce rapidly convergent approximations to the truncation errors of some slowly convergent in- 
finite series like the Dirichlet series for the Ricmann zeta function. Although the Euler-MacLaurin 
formula is extremely powerful in some cases, it is an analytic convergence acceleration method 
and cannot be used if only the numerical values of some elements of a slowly convergent sequence 
are known. Sequence transformations also try to accomplish an acceleration by constructing ap- 
proximations to the truncation errors of the elements of a slowly convergent sequence which are 
subsequently eliminated from the input data. However, they do this by by purely numerical 
means. 

Two practically very important classes of sequences are linearly and logarithmically conver- 
gent sequences. The partial sums of a power series with a nonzero, but finite radius of convergence 
constitute a typical example of a linearly convergent sequence, whereas the partial sums of the 
Dirichlet series for the Riemann zeta function, which is notorious for slow convergence, is a typ- 
ical example of a logarithmically convergent sequence. Thus, the properties of these two classes 
of slowly convergent sequences differ considerably, and so do the sequence transformations which 
are able to accelerate effectively either linear or logarithmic convergence. In general, the accel- 
eration of logarithmic convergence is more difficult than the acceleration of linear convergence, 
both theoretically and practically. 

Two powerful accelerators for linear convergence are Aitken's iterated A 2 process and Wynn's 
cpsilon algorithm which, in the case of the partial sums of a power series, produces Pade approx- 
imants. However, these two transformations fail to accelerate logarithmic convergence, whereas 
the powerful accelerators for logarithmic convergence - for example Richardson extrapolation or 
Wynn's rho algorithm and its iteration - fail to accelerate linear convergence. 

There are also some other problems with the acceleration of logarithmic convergence. In the 
literature on convergence acceleration, the transformations mentioned above are usually treated 
in their so-called standard forms. However, the standard forms are not able to accelerate the 
convergence of those logarithmically convergent sequences whose elements decay like a negative 
nonintegral power of the index. There is no a priori reason why sequences with a nonintegral 
decay parameter cannot occur in the context of oligomer calculations. Therefore, we also discuss 
transformations like Osada's variant of Wynn's rho algorithm or the weighted A 2 algorithm 
introduced by Bj0rstad, Dahlquist, and Grosse, which are both able to accelerate the convergence 
of sequences of that kind. 

We carried out calculations of the Hartree-Fock ground state energy for oligomers H — 
(CH=CH)n — H using a ST0-3G basis set with a fixed geometry, and extrapolated these en- 
ergies to the infinite chain limit corresponding to all-trans polyacetylene. The extrapolations 
were greatly facilitated by the fact that no accelerators were needed for the more troublesome 
logarithmic convergence. Contrary to previous belief we showed that the truncation error for 
the average energy per monomer unit contains linear term in 1/-/V, where N is the number of 
monomer units, but no higher powers. This leading contribution can be eliminated easily by 
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forming differences, and the higher order contributions converge exponentially. Thus, accel- 
erators for linear convergence sufficed for our purposes and Wynn's epsilon algorithm proved 
very effective in this regard. The ground state energy is just one property of interest. Future 
applications will deal with nonlinear optical properties and electronic transition energies which 
have a different dependence upon N and are more problematic from the practical point of view 
since they converge more slowly. 
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